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Abstract 

In this paper, we introduce a high-order accurate constrained transport type finite 
volume method to solve ideal magnetohydrodynamic equations on two-dimensional tri- 
angular meshes. A new divergence- free WENO-based reconstruction method is devel- 
oped to maintain exactly divergence-free evolution of the numerical magnetic field. A 
new weighted flux interpolation approach is also developed to compute the z-component 
of the electric field at vertices of grid cells. We also present numerical examples to 
demonstrate the accuracy and robustness of the proposed scheme. 



1 Introduction 

The ideal MHD equations model the dynamics of an electrically conducting fluid. Numerical 
solutions to magnetohydrodynamic (MHD) equations are of great importance to many ap- 
plications in astrophysics and engineering. Many efforts in solving the ideal MHD equations 
numerically have focused on the divergence-free evolution of the magnetic field implied by 
the induction equation 

— + VxE = 0. (1.1) 

Here B is the megnetic field, and E is the electric field defined by E = — u x B for ideal 
MHD. u is the velocity. J = V x B is the current density. The induction equation ensures 
that the magnetic field remains divergence-free if it is divergence-free initially. In numerical 
simulations, maintaining discrete divergence-free is also important. Previous studies [HI E] 
have shown that a divergence error on the order of numerical truncation error introduced by 
the numerical scheme can lead to spurious solutions and the production of negative pressures. 
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To name a few methods to ensure divergence-free evolution of the magnetic field, these in- 
clude Hodge projection approach f39], Powell's source term formulation [30], locally divergence- 
free discontinuous Galerkin (DG) method |25l E], constrained transport (CT) methods 
[38| IT2| I3H [T5| [32| El El |5l [19], generalized Lagrange Multiplier method [16], and many 
others [IIl|2l[37]. 

Despite these advances, almost all previous works have been focused on structured 
meshes. The CT type divergence-free formulation on structured meshes has been achieved 
at the second order accuracy in [H |5] and higher order accuracy in [8]. Several problems 
with complex geometry require the use of unstructured meshes. It is, therefore, desirable to 
design high order accurate divergence-free formulation for unstructured meshes. 

For the CT type formulation on structured meshes, the second-order accurate represen- 
tation of the magnetic field at the cell center can always be obtained by averaging the facial 
magnetic field. However, for the unstructured meshes, this is much harder to do, as there 
is no concept of arithmetic averaging of facial magnetic field to the center of the grid cells. 
As a result, the zone averaged magnetic field has always to be obtained via a reconstruction 
process on unstructured meshes. This makes divergence-free MHD on unstructured meshes 
slightly more challenging than the same process on structured meshes. 

In this paper, we introduce a divergence-free WENO reconstruction-based finite volume 
scheme up to the third order accuracy for solving ideal MHD equations on two-dimensional 
triangular meshes. ENO and WENO finite volume schemes have been introduced in many 
previous works for solving scalar conservation laws as well as compressible hydrodynamical 
flow problems using unstructured meshes [201 E ESI 1211 1221 [221 HI] • However, to the best of 
our knowledge, divergence-free high order (> 2) accurate finite volume schemes for solving 
ideal MHD equations on triangular meshes have not yet been available. To satisfy the 
divergence-free constraint of the magnetic field, we employ the CT framework. The basic 
idea of the CT framework adopted in the present paper is to introduce a staggered magnetic 
field at cell edges in two spatial dimensions (2D) (or faces in three spatial dimensions) and 
a staggered electric field at cell corners (or edges in three spatial dimensions) so that the 
computed magnetic field conserves a discrete definition of the divergence. To achieve this, 
a weighted flux interpolation approach based on [3] is introduced in this paper to compute 
the 2;-component of the electric field. To achieve high order accuracy, a new divergence- 
free WENO reconstruction method is introduced to reconstruct a cell centered magnetic 
field from the staggered allocated magnetic field on cell edges in two spatial dimensions. 
Additionally, the reconstructed piecewise smooth magnetic field is consistent at a cell edge 
by having the same cell edge-length-averaged value of normal component of the magnetic 
field when evaluated by using reconstructed magnetic field supported on triangles sharing 
this edge respectively. For the cell centered variables, the WENO reconstruction described 
in [221 HZ] is utilized. Numerical experiments show that the present divergence-free WENO 
reconstruction-based finite volume scheme is robust and accurate. 

The paper is organized as follows. Section [2] describes the CT type finite volume formu- 
lation to solve the ideal MHD equations. We start with introducing governing equations, 
notations for domain partition and discretization. Specifically, the proposed weighted flux 
interpolation approach to compute the 2;-component of the electric field is described in sub- 
section 2^ Section |3] describes the proposed reconstruction algorithm. The second-order 



accurate and the third-order accurate divergence-free WENO reconstruction methods are 
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catalogued in detail in subsection 3A_ Numerical tests are given in Section |4] to demonstrate 
the accuracy and non-oscillatory properties of the proposed scheme by computing smooth 
solution and shock wave related problems. We draw conclusions in Section |5] 



2 Finite Volume Formulation 

Ideal MHD governing equations in the conservation form can be expressed as 

dtV + a,F(U) + dyGiV) = , 



where 
and 
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Here p = Pgas + B ■ B/2 is the total pressure, 
following equation of state 



is the gas pressure that satisfies the 



(7-l)(e--pu-u--B-B) 



with u 



\Uxj Uy, Uz J 



and B = (B^, By, B^) . For a 2D ideal MHD problem, we have 

Ez = -UxBy + UyBx . (2.4) 

We employ the CT approach and the Godunov type finite volume scheme to solve Eq. 



(2.1). To this end, the physical domain is partitioned into a collection oi J\f triangular 
cells JCi so that f2 = IJi^i '^^ define 

rh = {}Ci:i = l,--- ,X} . (2.5) 

We also collect cell edges Cj to form 

(Eh = {Cy. J = !,■■■ ,m , (2.6) 

where TVg is the total number of edges in the partition. For every cell edge Cj, we uniquely 
identify an edge unit normal and tangent Cj- Here (^j is obtained by rotating rij 90 degrees 
in the counterclockwise direction. For simplicity, we assume that there are no hanging nodes 
in the partition Th- Let the edges of cell /Cj be denoted as dlCi^i, i = 1,2, 3. For convenience 
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in discussion, we define a mapping between the local cell edge index ^ of cell /Cj and the 
global edge index j such that 



i = (ii{3) and J=C'W- 
We also define the mesh parameter h to be 



h|c^ = the diameter of /Ci 

h = max,c.e77, ^x;. 



the longest side of /Cj 



(2.7) 



(2i 



We place the magnetic field variables -B^,. and By at the cell edges to maintain the global 
divergence-free evolution of the magnetic field; the z-component of the electric field at 
the cell vertices; and the conservative variables p, pu and e and on the cells. Br^ and 
By are always initialized to be divergence-free. The Godunov type finite volume scheme is 
utilized to evolve p, pu, e and Bz on the cells and the normal component of the magnetic field 
within the xy-plane on the cell edges. To evaluate at cell vertices, the flux-interpolated 
approach introduced by Balsara and Spicer ^ is further developed here. 

For convenience in discussion, we introduce notations U''^ = {p, pu,6, B^)^ and B^^ = 
{Bx,By) so that 

9iU^ + 9,F^(U) + 9,G^(U) = , 



(2.9) 



where 
F^(U) 

And 
where 
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(2.10) 



dtW^ + a,F^(u) + a,G^(u) = 



F^(U) 



G^(U) 



UxBy -\- UyBx 





(2.11) 



(2.12) 



Thus solving Eq. (2.1) is equivalent to solving equations (2.9) and (2.11) together. 



2.1 Semi-discrete finite volume scheme for the cell-centered U 



H 



Taking the cell /Cj 



, A/", in partition (2.5) as a discrete control volume, the semi- 



discrete finite volume method for solving Eq. (2.9) is formulated by integrating (2.9) over 
the cell /C, : 



Ufc,(t) + 



dt 



Ff,Gn -11,^1 = 



(2.13) 



dKi 
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where Vj^^it) is the cell average of the k^^ {k = 1, ■ ■ ■ , 7) component of on /Cj, is the 
f^th component of F'^, Gjf is the k^^ component of G^, and rij is the outward unit normal 
of the boundary of the cell /Cj. |/Cj| is a shorthand notation for the area of /Cj. 

To solve Eq. (2.13) numerically, we evaluate the flux integral by Gaussian quadrature 
rule with the exact value of (F^, G^) -rij being replaced by the Lax-Friedrichs flux F*{x, y, t) 
given by 

ni^,y,t) = I [(Ff(U-),Gf(U-)) + (Ff(U+),Gf(U+))] .n,-|(Uf+-Uf-). (2.14) 

Here a is taken as an upper bound for the eigenvalues of the Jacobian in the rij direction; U~ 
(or U^'~) and (or U'^'"'") are the numerical values of U (or U^) inside the triangle and 
outside the triangle at the Gaussian point. To this end, we obtain the following semi-discrete 



finite volume scheme for solving Eq. (2.9) 



■Ti^lkM + W-\ F*,cir = 0, (2.15) 
J die. 



where U;^ j(t) is the approximate cell average of the k component of U on the cell /Cj. 

2.2 Semi-discrete finite volume scheme for the edge-centered nor- 
mal component of B^^ 

The 2D constrained transport scheme developed in the present paper is based upon cell 
edge-length-averaged magnetic field located at the edges of grid cells. On every cell edge 



Cj G we solve Eq. (2.11) to evolve the normal component of B^^ with respect to the 
defined cell edge unit normal . Denote the normal and tangential contribution of B^^ in 
directions given by rij and C,j to be Bn and -B^ respectively. We rewrite Eq. (2.11) by Bn 
and B(^ to obtain 

Here m„ and are the components of velocity u in the and C,j directions respectively. 
Let Bnj be the edge-length-averaged Bn on the edge Cj defined by 

/" 5nrfC, (2.17) 



where is a shorthand notation for the length of the edge Cj. Integrating Eq. (2.16) 



along the cell edge £j, the semi-discrete finite volume scheme to evolve Bn,j numerically on 
Cj can be expressed as 

^5;— = -^^^^4^^^^, (2.18) 

dt ' \Cj\ ' ^ ' 



smce 



= -UnBc_ + Uc^Bn ■ 
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Here -B/i,nj is the approximate cell edge-lengtli-averaged -B„ on Cj. Ez{Cj^e) is numerical 
approximation of the z-component of E at the end point of Cj, and Ez{Cj^s) is numerical 
approximation of the 2;-component of E at the starting point of Cj. In the direction of ^j, 
the two end points of the edge Cj are defined to be the starting and the end point of Cj 



respectively. The method to compute is described in Section 2.3 



2.3 Computing Ez at the vertices of cells by flux interpolation 




Figure 1: The stencil to compute Ez at the vertex V. Cell edges Cq, - ■ ■ , £4 which separate 
triangular cells JCq, - ■ ■ , /C4 all have one end at V. no, ■ ■ ■ , n4 are unit normals of these edges 
respectively. 



In our scheme, one has to obtain the electric field E^ at vertices of the triangular mesh 
(see Fig. [T]). In [3], it was shown that there is a dualism between the electric field and 
the properly upwinded flux. In fluid dynamics, such a flux takes on contributions that 
are upwinded normal to a zone face. For MHD, the electric field at the vertex V in Fig. [T] 
should take on properly upwinded contributions from all possible directions. This necessarily 
would require a mult i- dimensional Riemann solver. For structured meshes, such a multi- 
dimensional Riemann solver has been presented in [9]. Unfortunately, a multi-dimensional 
Riemann solver that works for MHD on unstructured meshes has not been presented in the 
literature. For that reason, we use the available ideas on multi-dimensional upwinding from 
[3] and the idea of doubling dissipation in each direction from |27l HH] . 

Below we describe an algorithm to compute the z-component of the electrical field Ez 
at vertices of the mesh. The algorithm results in an upwinded choice of Ez in a multi- 
dimensional fashion. 

See Fig. [1} Suppose triangles /Cq, ■ ■ ■ ,/C4 meet at the vertex V. The edges shared by 
triangles are labeled by " " " ? -^4 and the associated unit normals of edges by hq, ■ ■ ■ , n4 
respectively. 



On the each edge Ci,l = 0, - ■ ■ ,4: = ny, using Eq. (2.4), which shows the dualism between 
Ez and flux, we obtain 

^2,/(a;v,Z/v) = -Flf,6 (U"(xv,i/v),U+(xv,z/v)) , (2.19) 
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from numerical flux interpolation. Here (xv,l/v) are the coordinates of V. F^^ is the Lax- 
Friedrichs flux with double dissipation for solving Eq. (2.1); and F^pfi is the 6^^ component 
of Fii;'. Let /CJ"* denote the interior of cell /Cj. 



int 



U {xv,yv) = lim(^,y)^(^^,j^^)U(x,2/,t) , where {x,y) e /C; 
lJ+{xv,yv) = lim(^x,y)^{xy,yy)U{x,y,t) , where {x,y) G . 

Thus 

F^^(U-, U+) = ^ [(F(U-), G(U-)) + (F(U+), G(U+))] ■ n, - a(U+ - U") . (2.20) 

Here a is taken as an upper bound for the eigenvalues of the Jacobian in the direction. 
If the flow is locally smooth, we can take the arithmetic average 

^ riv 

Ez{xv,yv) = ^ y]^2,«(a;v,2/v) 

^ 1=0 

to obtain a unique Ez at the vertex V. However, when discontinuities are present it is bene- 
flcial to allow the evaluation of to locally adjust to those discontinuities. To achieve this, 
we design switches to detect strong magnetosonic shocks and strongly compressive motions 
and the direction of propagation of the discontinuity. 

For this purpose, we flrst use a least square approach to construct linear proflles of 
pressure and velocity at the vertex V respectively as follows. See Fig. [TJ Briefly, we flrst 
compute on ICi, / = 0, ■ ■ ■ ,4, pressure Pgas,i and velocity {ux,i,Uy^i) from cell average values 
of conservative variables at the cell centers. Let v = {pgas,Ux,Uy)'^', and the s*'^ component 
Vs{x, y), s = 1,2, 3, of V be represented by a linear polynomial 

Vs{x, y) = vo,s + —J^i^ ~ ^v) + —f^iy ~ yv) (2-21) 

where h is the mesh parameter; {A^Vs, AyVg) is the undivided difference approximation to 
the gradient of the exact proflle of pressure or velocity at V. Parameter values Vq^s, ^x^s, 
AyVs are determined by solving 

^xVs AyVs 

vo,s H — ^(^' - xv) H — j^yyi " y^i = ^^.^ ' / = o, ■ ■ ■ , 4 

in the least square sense. Here {xi,yi) are the coordinates of the cell center of /C;; Vs^i stands 
for the value of the s*'* component of v on cell K-i, which are Pgas,i, Ux,i and Uy^i respectively. 

The flrst switch, SWl, which is used to pick out strong magnetosonic shocks or conflgu- 
ration that may develop into such a shock, is accomplished by taking the undivided gradient 
of the pressure at the vertex V and comparing it with the minimum pressure in the vicinity. 
SWl is switched on if 

\A^Pgas\{Xv,yv) + \^yPgas\{Xv,yv) > (3 iam{pgas,0, " " " ,Pgas,nJ {2.22) 

and is switched off otherwise. Here we use (3 = 0.5. \Ax\{x\;,y\;) and | Aj,|(xv, l/v) are 
shorthand notations for absolute values of A^; and Aj^ at (xv,|/v) respectively. 
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The second switch, SW2, which is used to pick out strong compressive motions at the 
vicinity of the vertex V, is accomphshed by comparing the undivided divergence of the 
velocity to the smallest local signal speed. SW2 is switched on if 

- 5min(Co, ■ ■ ■ , C„ J > (A^..^^. + AyUy){xv,yv) (2.23) 

and is switched off otherwise. Here we use 6 = 0.1. 
On cell /C/, 

Ci = \ , / = 0, ■ ■ ■ 

\ Pi Pi J 

When either SWl or SW2 is switched on, it means that the region has a shock in it. In 
this case, we need to pick out the direction along which we want to upwind the evaluation 
of the electric field. We use a weighted combination to do that. 

We estimate the direction ns = {ns,x,ns^y)'^ of the strong shock in the vicinity of the 
vertex V by 

flg^ = ^xPgas ^2 24) 

Then for each Ez^u we compute the associated weight wi by 

= ^^S^ (2.26) 

where ai is defined by 

ai = (ns • niY + IQ-^ . (2.27) 



Here the small number 10 ^ is to avoid division by zero in Eq. (2.26). 
We obtain 

Ezixv,yv) = ^ujiEz,i{xv,yv) (2.28) 

1=0 

at vertex V when discontinuity or strong compression is present. 
2.4 Time discretizations 

The method of lines approach is used to evolve the solution on the triangulated domain. 
Specifically, the third-order accurate TVD Runge-Kutta method [S] is used to solve ordinary 



differential equations (2.15) and (2.18) 



3 WENO-based Reconstruction 

The main ingredient of a high order accurate finite volume scheme is a reconstruction algo- 
rithm, which reconstructs a smooth and high degree polynomial approximation of solutions 
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from average values computed by the base finite volume scheme at the end of every Runge- 
Kutta stage. In return, the reconstructed polynomial is used for evaluating numerical fluxes 
in the subsequent calculation. In this section, we solve the following two sub-problems of 
reconstruction: 



Sub-problem 1. Given edge-length-averaged normal component B^j of B^^ define on cell 
edges Cj G and a positive integer g, for each cell /Cj, reconstruct an essentially non- 
oscillatory and divergence-free magnetic field B^^ supported on /Cj. Here Bf^ G Pg(/Ci)^. 
PqiJCi) is the space of polynomials of degree at most q supported on /Cj. B^^ is a (g + 1)*'^ 
order accurate approximation to exact B^^ (when it is smooth) on cell /Cj. Moreover, 



I , Bf^ • n,,dr = , £=1,2,3. (3.1) 



Here dlCi^ is the t edge of cell /Cj; rij ^ is the associated unit normal of this edge; \d]CiA 
is its length. B^^-ij-^^ is the edge-length-averaged normal component of the magnetic field 
on the cell edge (9/C^. The local edge index i a nd th e global edge index j is related by the 



mapping function (2.7). We note that condition (3.1) implies that the piecewise smooth B^^ 
agrees at the adjacent cell edges by edge-length-averaged mean values. For our proposed 
scheme, we apply the reconstruction algorithm to solve this sub-problem at the end of every 
Runge-Kutta stage by using mean values Bh,n,j, which is the numerical approximation of 

Bn,j ■ 

Sub-problem 2. Given cell average values Vi of a function v{x,y) on each cell /Cj and 
a positive integer q, for each cell /Cj, reconstruct an essentially non-oscillatory polynomial 
Pi{x,y) of degree at most q which has the mean value Vi and is a (g + l)*'^ order accurate 
approximation to v{x,y) on /Cj (when v{x,y) is smooth). For our ideal MHD problem, at 
the end of every Runge-Kutta stage, this problem is solved for Vi replaced by the cell average 
values of the every component of computed by the base finite volume scheme. 

Before we describe the algorithm to solve these two reconstruction problems, we first 
recall several relevant concepts which will be used later in this section. See also [22] for 
details of related discussion. The level-0 von Neumann neighborhood of a triangle JC E Th 
contains the edge adjacent neighbors of /C and is defined to be the set 

OT°(/C) = 1^ e {/C} : ^ n /C is an edge of /c} . 

Here we neglect the subscript "z" of cells for convenience. The level-r von Neumann neigh- 
borhood is defined by the recursive definition 

<yr(IC) = j IJ ^'-^jC) I \ {/C} , for r > 1 . 

For instance, the level- 1 von Neumann neighborhood of /C is by merging level-0 von Neumann 
neighborhoods of cells in ^^"(/C) which are edge adjacent neighbors of /C. 

A critical component for the success of reconstruction on triangular meshes is the selection 
of a set of admissible stencils. Generally, this set should contain isotropic (or centered) 
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stencil for achieving good approximation in smooth regions, and anisotropic (or one-sided 
and reverse-sided) stencils to avoid interpolation across discontinuities [221 [T7] . 

In order to construct these anisotropic stencils, the sector search algorithm [221 [20] is 
utilized to construct forward sectors as well as backward sectors. Fig. 2(b) shows three 
forward sectors FSg, s = 1, 2, 3 of cell /Cq. A forward sector is spanned by a pair of edges 
of /Co- 



Fig. 2(c 



shows three backward sectors BSs,s = 1,2,3 of JCq. A backward sector 
is defined by having its origin at the midpoint of a edge of JCq and its two boundary edges 
passing through the other two midpoints of remaining edges of /Co- 

When we perform a WENO reconstruction for a cell /C G 7/i, in addition to construct a 
central stencil, within each of the sectors of /C, we construct either an one-sided stencil (when 
it is a forward sector) or a reverse-sided stencil (when it is a backward sector). Additionally, 
when we construct an anisotropic stencil, we only include cells in von Neumann neighbors 
of /C, whose barycenters lie in the corresponding forward or backward sector. 



3.1 Divergence-free WENO-based reconstruction for B^^ on cell 

Here we describe a new divergence-free WENO-based reconstruction strategy for the mag- 
netic field on triangular grids based on our recent work [10] and [HI [7]. This solves the 
Sub-problem 1. We require that the reconstructed B^^ supported on cell /C, must satisfy 
the divergence-free condition on /Cj internally, is a (g -|- 1)*'' order accurate approximation to 
exact B"^^ on /Cj and also retains consistency at the cell boundaries in the sense defined by 



Eq. (3.1). We summarize the reconstruction algorithm as follows: 



Step 1. For every grid cell /Cj, we identify a set of admissible reconstruction stencils = 
{T^™-' : m = 1, ■ ■ ■ ,7} using the method introduced in [201 1221 [IZ]- Here T^^^ is the 
central stencil; xj^^ and Tj^^ are the one-sided stencils constructed in the forward 
sectors of /Cj; and T^^^ and T^^^ are the reverse-sided stencils constructed in 

the backward sectors of /Cj respectively. This choice of stencils allows to better limit 
oscillations of polynomial approximation of solutions supported on /Cj. 

Step 2. We then use every stencil to reconstruct preliminarily a divergence-free magnetic field 
Bf'("^\m = I,-- - ,7, with every component of B^^'^™''' represented by a polynomial 
function from cell edge-length-averaged values of normal component of the magnetic 
field B^^ defined on edges of cells contained in the stencil. 

Step 3. For each preliminarily reconstructed B^^'*'™'', a smoothness indicator Um is computed 
with Ylln=i^m ~ The final nonlinearly stabilized WENO reconstruction B^^ is 
defined by a weighted combination of Ylln=i ^m^i^'^^^ and is exactly divergence- free. 



3.1.1 The second order accurate reconstruction for 

To explain this reconstruction algorithm clearly, we first describe the second order accurate 
reconstruction algorithm. We note that our reconstructed B^^ belongs to {Pi(/Ci)^, V-B^^ = 
0}. The following linear polynomial expression is employed for representing the preliminarily 
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Figure 2: Stencils for reconstructing second order accurate cell-centered divergence-free mag- 
netic field on cell /Cq. Normal components of the magnetic field on solid line edges are utilized, 
(a) The central stencil, (b) Three forward sectors FSi, FS2 an FS3 of cell /Cq formed by 
spanning a pair of edges of /Co respectively. The cells of three one-sided stencils formed in 
each of the forward sector is shown here, (c) Three backward sectors BSi, BS2 and BS^ of 
/Cq. a backward sector is defined by having its origin at the midpoint of an edge of /Co and 
its two boundary edges passing through the other two midpoints of remaining edges. The 
cells of three reverse-sided stencil formed in each of the backward sector is shown here. Note 
that the same notations are used for different normal components of the magnetic field and 
cells in (a), (b) and (c) to avoid introducing too many notations. 

reconstructed B^^'-M = (B^r\ ^i""^)^ as weU as Bf 

-Si"*^ (X, y) = ao^rn + ai^rnX + a2,r„l/ , 
B^r'' ix,y) = + hi^rnX + b2,my ■ 



(3.2) 



2| shows stencils used to reconstruct Bq^ on cell JCq. The central stencil is shown 

= 1,2,3 of /Co as well as cells 
shows three backward sectors 

The 



Here we drop the subscript "i" of cells for convenience. The divergence- free condition V ■ 
■Qxy,{m) _ g gjygg ^]-^g equatiou 

ai,m + &2,m = , (3.3) 

for linear polynomial representation of B^^ by matching coefficients. This reduces one of the 
degrees of freedom in preliminarily reconstructing B^^'*^™^ 
Fig. 

in Fig. 2(a) Fig. 2(b) shows three forward sectors FSs,s 
used to construct corresponding one-sided stencils. Fig. 2(c 

BSs, s = 1, 2, 3 of /Co and cells chosen to construct corresponding reverse-sided stencils 
details of constructing one-sided and reverse-sided stencils in these sectors and reconstructing 
linear polynomial Bq^ supported on /Cq are given below. 

Computation of the second order reconstruction polynomial for central stencil. 

See Fig. 2(a) The central stencil consists of /Cq, and its three adjacent neighbors /C^, s = 
1, 2, 3, which are in the level-0 von Neumann neighbor of /Cq. 

On /Co 5 we arbitrarily choose two cell edge- length- averaged values of normal components 
of B^y, denoted by Bq and Bi. On each of /C^, s = 1,2,3, we choose one cell edge-length- 
averaged value of normal components of defined on the cell edge which is not the common 
edge between /Co and /C^, and denote these average values by B2, B3 and B4 respectively. The 
cell edges on which these values are defined are relabeled by Cs '■ s = 0, ■ ■ ■ ,4 for convenience. 



The solid lines in Fig. |2 (a) [ indicate these edges. We then solve 



ai,r 
1 



'2,m 



J^^ ^xy,im) . = Bs, S = 0, 



(3.4) 



to obtain a candidate B^^'(^\ Here m = 0; \Cs\ is the length of Cg] and n^ is the unit normal 
of Cs- 

Notice that our choice of cell edges is such that we never form closed loops of edges. This 
is because that the divergence-free condition ensures that each closed loop has one redundant 



piece of information. Also notice that Eq. (3.2) only has five independent degrees of freedom 



and we have selected five edges which do not close from the central stencil. We remark that 
this general principle also applies to reconstructing B^^ on other stencils. 

Finally we point out that for this central stencil, there are other possible ways to choose 
cell edges for reconstruction. See Fig. 



2 a 



For instance, we could use values defined on the 
dashed line edges from every neighboring cells of /Co together with ones defined on the solid 
line edges of /Co to reconstruct the polynomial. This is also the case for reconstruction by 
using other stencils. However, we notice that as long as we follow the principle of choosing 
edges close to /Co, the results are not sensitive to choices of edges. 

Computation of the second order reconstruction polynomials for one-sided sten- 
cils. 



Fig. 2(b) shows forward sectors and corresponding one-sided stencils utilized for recon- 
struction. In Fig. 2(b), cells /Ci, /C2 and /C3 are three neighbors of /Co; and JCso, /C^i are 



two neighbors of ICg (other than /Co), s 



(2) 

1,2,3. We form three one-sided stencils: Tg 
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{/Co,/C3,/C3o,/C3i} in FSi, T^^ = {/Co, /Ci, /Cio, /Cn} in FS2] and T^' = {/Co, /C2, /C20, /C21} 



(4) 



in 



With every one-sided stencil T 



(m) 
B ■ 



m = 2,3,4, we will utilize cell edge-length-averaged 
values of normal component of the magnetic field B^ ^ on cell edges to reconstruct prelim- 



inarily a divergence-free polynomial in the form of (3.2). Take stencil T^^ for example. 



On /Co, we utilize two average values of normal component of B^^ defined on its two edges 
respectively, and denote them by Bq and Bi. On /C3, we use one average value of normal 
component of B^^ defined on the cell edge which is not the common edge between /Cq and 
/C3, and denote it by i?2- On each of /C30 and /C31, one average value defined on the cell edge 
which is not shared by /Ci and /C30 (or by /C3 and /C31) is employed respectively, and denote 
them by B3 and B4. 

We now use {Bg : s = 0, ■ ■ ■ , 4} to reconstruct preliminarily a piecewise linear B"^^'*^^^ by 
solving Eq. (3.4). We similarly compute on the other two one-sided stencils to obtain two 



candidates respectively, denoted by B^'^'(^\ and B^'^'^'^^ 



Computation of the second order reconstruction polynomials for reverse-sided 
stencils. 

The reverse-sided stencils are constructed by using cells within the backward sectors. 



Fig. 2(c) shows reverse-sided stencils which are constructed in the backward sector BSg, s = 
1,2,3 for reconstructing on /Cq- Here cells /Ci, /C2 and /C3 are three neighbors of /Co; 
and ICso, ^si are two neighbors of fCs (other than /Co), s = 1,2,3. We construct three 



reverse-sided stencils: T^'' = {/Co, /Cio, /C20} in BSi, T^' = {/Co, /Cn, /C31} in -85*2; and 



(6) 



T^^ = {/Co,/C2i,/C3o} in 55*3. 

Unlike the central or one-sided stencil case, here we employ a constrained least square 
method to solve the reconstruction problem by using each of the reverse-sided stencils. Due 
to the divergence-free condition, in principle, we only need 5 additional conditions to uniquely 



determine functions (3.2). However, a reverse-sided stencil can provide 6 admissible average 
values of normal component of B^^ (or 6 conditions to determine (3.2 )). See Figure 2(c) Take 



the stencil T^^"* for example. The edges on which the defined average values are employed for 
preliminary reconstruction are indicated by solid lines. On /Cq, we can use two edge values 
Bq and Bi. On /Ciq, we can use B2 and B^, and on /C205 we have i?4 and B^ to use. To 
avoid a bias in choosing values from we use all of these 6 average values and solve the 
following constrained least square problem for the preliminary reconstruction: 



\C4 



B. 



(3.5) 



subject to: 



_L 
\Ci\ 



+ h 



2,m 



_ J^^ Bxy,{m) . ^^^^ 



0,1 



(3.6) 



When solving this constrained least square problem, Eq. (3.6) is satisfied exactly; while 
Eq. (3.5) is satisfied in the least square sense. To improve the divergence-free aspect of 



the solution to this constrained least square problem on stencil cells other than /Cq, one 
can substitute Eq. (3.3) into Eq. (3.5) before solving (3.5)-(3.6). This ensures that the 
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preliminarily reconstructed magnetic field on the cell of interest, JCq, is divergence-free, and 
matches the magnetic field defined on the cell edges bounding /Co by mean values exactly. 
Additionally, the preliminarily reconstructed magnetic field is the best approximation to the 
magnetic field on other cells in the reverse-sided stencil. 



We denote candidates obtained by solving equations (3.6)-(3.5) in a constrained least 
square manner for every reverse-sided stencils B^^'^^\ 'Q^vA^) and B^'^'^^^^ respectively. 

Computation of weights for second order WENO reconstruction. 

We now apply a weighted combination of {B^^'^™) : m = 1, ■ ■ ■ , 7} to finalize Bq^ using 



the idea of WENO [351 HE]. Let B^^'-^'") be expressed by Eq. Q. 

We compute a quantity a^, which is the reciprocal of a smoothness measure by 



am. = : : r, ri^ r, tit , m = 1, ■ ■ ■ , 7. 



Here we take e = 10~^ to avoid division by zero. 6^ = 10 when m = 1; and 6^ = 1 otherwise. 
This follows the idea in [23]. The parameter bm allows to have more weight on the central 
stencil, which provides better accuracy when the solution is smooth. 
The weight Um is computed by 



Finally, we reconstruct the piecewise linear polynomial approximation Bq^ on /Cg by 

7 

B^ = J]]a;„B^^'(™). 

m=l 

Notice that Bq^ satisfies the divergence-free condition: V ■ Bq^ = exactly. This completes 
the second order accurate reconstruction for approximating B^^ on /Cq- In our proposed 
scheme, at the end of every Runge-Kutta stage, we apply this reconstruction strategy by 



using values -B/i,n,j computed by the base finite volume scheme (2.18 ) to do the reconstruction 



3.1.2 The third order accurate reconstruction for B^^ 

The third order accurate divergence-free reconstruction of B^'^ on cells is a straight forward 
extension of the second order accurate reconstruction described in Sec. I3.1.1[ Thus we 
catalogue the key steps in implementing the third order accurate case in this subsection. 

The reconstructed Bf^ supported on cell /Cj, fCi G Th, belongs to {P2(/Cj)^? ^ ' = 0}. 
To avoid introducing too many notations, representation of the preliminarily reconstructed 

•Qxy,(m) _ (j^^^^^ By"^^) and B^^ is redefined by the following quadratic polynomial expres- 



sion 

Bi"'\x, y) = ao,m + ai^rnX + a2,my + a^^mx"^ + a^^mxy + a^,my'^ , ^2 

Bj^^ {X,y) = 6o,m + bi^rnX + b2,my + bs^rnX'^ + b^^^Xy + 65,m2/^ ■ 
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(c) 

Figure 3: Stencils for reconstructing the third order accurate cell-centered divergence-free 
magnetic field as well as for reconstructing the third order accurate polynomial approxima- 
tions for cell centered variables on cell JCq. To reconstruct the divergence-free magnetic field, 
normal components of the magnetic field on solid line edges are utilized, (a) The central 
stencil, (b) Three forward sectors FSi, FS2 an FS3 of cell /Cq formed by spanning a pair of 
edges of ICq respectively The cells of three one-sided stencils formed in each of the forward 
sectors are shown here, (c) Three backward sectors BSi, BS2 and BS3 of /Cq. A backward 
sector is defined by having its origin at the midpoint of an edge of /Cq and its two boundary 
edges passing through the other two midpoints of remaining edges of /Co- The cells of three 
reverse-sided stencils formed in each of the backward sectors are shown here. Note that the 
same notations are used for different normal components of the magnetic field and cells in 
(a), (b) and (c) to avoid introducing too many notations. 

The subscript "i" of cells is also dropped for convenience. We note that the divergence-free 
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condition V ■ B^f'(™) 



gives the following three equations 



,m 








(3.8) 



for quadratic polynomial representation of B^'^ by matching coefficients. This reduces three 
of the degrees of freedom in preliminarily reconstructing the third order accurate B^y'(^\ 

Fig. |3] shows stencils used for the third order accurate divergence-free reconstruction. The 
central stencil is shown in Fig. 3(a) Fig. 3(b) shows three forward sectors FSs,s = 1,2,3 
and cells in these forward sectors used to construct one-sided stencils. For the convenience, 
we still denot e the one-sided stencil in BSi by Tj^\ the one in BS2 by T^^^ and the one in BS3 

1,2,3 and cells in these backward 



Fig. 



3 c 



shows three backward sectors BSg, s 



sectors used to construct corresponding reverse-sided stencils. Similarly, the reverse-sided 
stencil in BSi is denote by T^^\ the one in BS2 by T^^^ and the one in BS3 by T^^^ The 
details of constructing these stencils and reconstructing polynomial B^^ are given below. 



Computation of the third order reconstruction polynomial for central stencil 

See Fig 



3(a, 



The central stencil T^^ of /Cq contains cells in level-1 von Neumann 
neighbors of /Co and /Co itself. On cell /Cq, we arbitrarily choose two cell edge- length- averaged 
values of normal components of B^^. On each of the remaining cells, we choose one edge- 
length- averaged values on the edge which is connected to one of vertices of /Cq. We relabel 
these average values by {B^ : s = 0, ■ ■ ■ , 10}. We then solve the following constrained least 
square problem to obtain B^^''^^-' 



{ jtl k B^^'^"^^ ■ ^sdr = Bs, s = 2, 



10 



(3.9) 



subject to: 



0'l,m + &2,»n 
2a3,m + ^4,m 
04 + 265 

1 r B^'2^'("*) ■ nidr 



, 
, 


Bi, 



(3.10) 



/ = 0,1 . 



Here m = 1. We remark that when solving this constrained least square problem, Eq. (3.10) 
is satisfied exactly; while Eq. (3.9) is satisfied in the least square sense. Similar to the 



second order accurate preliminary reconstruction on reverse-sided stencil, to improve the 
divergence-free aspect of the solution to this constrained least square problem on stencil 
cells other than /Cq, Eq. (3.8) is substituted into Eq. (3.9) before solving (3.9)-(3.10). This 



ensures that the preliminarily reconstructed magnetic field B"^'^'*^^) on /Co is divergence- free, 
and matches the mean values of the magnetic field defined on the cell edges enclosing /Cq 
exactly. 



Computation of the third order reconstruction polynomials for one-sided stencils. 
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Within each of the forward sectors, we construct an one-sided stencil. See Fig. 3(b 



Take the one-sided stencil T^^ constructed in the forward sector FSi for example. T^^ 
contains cells in level-3 von Neumann neighbors of ]Cq, which are close to the edges of -F^i, 
and /Co itself. Thus T^^^ = {/Co, /C3, /C30, /C31, /C4, /C5, /Ce, One-sided stencils T^^ and 
T^"^ are constructed similarly. 

To reconstruct quadratic polynomial B^^'*^^) on /Co by using Tj^\ we choose two cell edge- 
length-averaged values of normal components of B^^ with one of them defined on the edge 
shared by cells /Cq and /C3 (level-0 von Neumann neighbor). Then we choose 7 edge values 



from the remaining cells. The solid line edges shown in Fig. |3(b) give one admissible selection 
of these 7 edge values. We relabel these edge-length-averaged values by {Bg : s = 0, • • • ,8}. 
Then the following linear problem is solved to reconstruct preliminarily a candidate B^^'*^^^ 



2a3,m + &4,m = , 

_1_ J^^ ^xy,im) . = B„ S = 



{3.11] 



Here m = 2. The other two preliminarily reconstructed polynomials B^^'^^^ by using T^^ 
and B^?''^^) by using T^^^ are computed in the same manner. 

Computation of the third order reconstruction polynomials for reverse-sided 
stencils. 



See Fig. 3(c) The reverse-sided stencils are constructed in the backward sectors re- 
spectively. Take the reverse-sided stencil T^^^ constructed in the backward sector BSi for 

(5) 

example. Tg contains two level-1 von Neumann neighbors /C41 and /C40, one level-2 von 
Neumann neighbor /C4 which is adjacent to both /C41 and /C40, one level-3 von Neumann 
neighbor /C5 which is adjacent to /C4 and two level-4 von Neumann neighbors /C50 and /C51 
which are adjacent to /C5 (other than /C4). The other two reverse-sided stencils T^^ and T^^ 
are constructed in the same manner. 

To reconstruct the quadratic polynomial B^^'*^^^ supported on /Cq by using stencil tI^\ 
we note that there are also multiple approaches to select cell edge-length-averaged values of 
normal component of B^^. Fig. 3(c)| shows one choice. Normal components of B"^'^ defined 



on the solid line edges in stencil Tg are utilized. We choose two cell edge- length- averaged 
values of normal components of B^^ on two edges of /Co whose comment endpoint is in 
the backward sector BSi. We also choose four cell edge-length-averaged values of normal 
component of B^^ from edges of /C40 and /C41, which are level-1 von Neumann neighbors of 
/Co. We then choose three cell edge- length- averaged values of normal component of B^^ from 
edges of remaining cells in the stencil. 

We relabel these edge-length-averaged values by {Bs : s = 0, ■ ■ ■ , 8}. We then solve the 



linear problem (3.11 ) to obtain B"^'^'*^^^ The other two preliminarily reconstructed polynomials 



B^y'(^) by using Tg and B^^'*^"^^ by using Tg are computed in the same manner. 
Computation of weights for third order WENO reconstruction. 
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We now apply a weighted combination of {B^^'*^™) : m = 1, ■ ■ ■ , 7} to finalize Bq^ using 
the idea of WENO [351 [IB]. Let B^f'(™) = (^Bi'^\x,y), B^""\x,y)^ be expressed by Eq. 



(3.7). 



We first compute smoothness measures of x-component and ^/-component of B^^'^'"-' re- 
spectively by 

Si(B^"^))= lj2 J h-''{D^B^"'\x,y)fdxdy\ . (3.12) 



Here B stands for either Bi"^\x, y) or By"^\x, y). 
The am is redefined by 



(m)> 
^2/ > 



(e + 5/(5,,('")) + 5/(5, 



e = 10 ^ is used to avoid division by zero, bm = 10 when m = 1; and b^ = 1 otherwise [23] . 
The weight Um for the third order accurate divergence-free reconstruction is computed 

by 



Finally, we reconstruct the piecewise quadratic polynomial approximation Bq^ on /Cq by 

7 



m=l 



Notice that Bq^ satisfies the divergence-free condition V ■ Bq^ = exactly. This completes 
the third order accurate divergence-free reconstruction for the magnetic field on /Co- Again, 
at the end of every Runge-Kutta stage, we apply this reconstruction method by using values 



Bh,n,3 computed by the base finite volume scheme (2.18) to do the reconstruction 



3.2 WENO finite volume reconstruction for 

Here we describe an algorithm to reconstruct polynomials of degree g = 3 from given cell 
averages to solve the Sub-problem 2. 

Let [xi, yi) be the coordinates of the barycenter of cell /Cj. We use the following monomial 
expression of a second degree polynomial Pi(x, y) supported on /Cj: 

Pi{x,y) = ao^i+ai^i{x-Xi)+a2,iiy-yi)+a3^i{x-x,i)^+ai^i{x-Xi){y-yi)+a5^i{y-yiY . (3.13) 

To reconstruct a polynomial function approximation to a function v{x,y) on cell /Cj from 
cell average values Vi of v{x, y), we also follow the reconstruction algorithm described in Sec. 



3.1 except that in Step 2 of the algorithm, we use cell average values here for solving the 
Sub-problem 2; and we do not require the reconstructed Pi{x,y) to be divergence-free. 

For the self-completeness of the paper, we briefly describe the WENO reconstruction of 
the second degree polynomial Po{x, y) on cell /Co- We refer to [I7l[T8l|22] for description of the 
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first degree polynomial reconstruction. See Fig. |3] for all stencils used in the second degree 
polynomial reconstruction. On each of the stencil T^"^\m = 1, ■ ■ ■ ,7, we first reconstruct 
preliminarily a polynomial Pq™'' {x, y) supported on /Co by solving a system of linear equations 
(or a constrained least square problem) respectively. 

The central stencil T*^^-* to reconstruct preliminarily a second degree polynomial P^\x, y) 



on cell /Co is shown in Fig. 3(a) T*^^) consists of cell /Cq, its three neighbors /Ci, /C2 and 
/C3, and /Cso and /C^i which are two neighbors of /C^ (other than /Co), s = 1,2,3. Thus 
T^^^ = {/Co, /Cs, /Cso, /Csi, s = 1,2,3}, which consists of the level-1 von Neumann neighbors 
of /Co and /Co itself. The coefficients of Pq^^ (x, y) are determined by solving the following 
constrained linear problem 

I Jf^^^PQ^\x,y)dxdy = \JCsi\vsi , s = 1,2,3; / = 0, 1; 

subject to: (3-14) 



T\'^T J ^ -^7 ^7 ^1 



where |/C| is the area of cell /C; Vgi is the cell average value defined on cell Kgi] and Vr is the 
cell average value defined on /C^. 

Within in three forward sectors FS^, s = 1, 2, 3, we construct three one-sided stencils T^'^\ 



T(3) and T(4) respectively. See Fig. |3(bj] In FSu T^^) = {/Co, i^s, ^30, /C31, /C4, /C5}. Thus 
T*^^) consists of /Cq, level-1 von Neumann neighbors of /Co in this sector, and two additional 
cells /C4 and /C5 which are neighbors of level-1 von Neumann neighbors in this sector (other 
than /C3). The other two one-sided stencils T*^^^ and T^^^ are constructed similarly. 

Then by using every stencil T^"^\m = 2,3,4, we reconstruct preliminarily polynomials 
Pl^\x,y) respectively by solving the following linear system 

/ Pt\x,y)dxdy=\^s{T^'^^)\vs , s = l,---,6, (3.15) 

7a3(t(")) 

where A,(T('")) e T^") is a cell in T'^'^\ m = 2,3,4; Vs is the cell average defined on 
A,(T('")), and |A,(T(™))| is the cell area of A,(T(™)). 

To further improve the robustness of the scheme, the reverse-sided stencils are also in- 



cluded. See Fig. 3(c), Within three backward sectors BSs,s = 1,2,3, we construct three 
reverse-sided stencils T^^\ T^^^ and T*^^^ respectively. 

In BSi, T^^^ = {/Co, /C4, /C40, /C41, /C5, /C50, /C51}. The method to construct this stencil 
is the same as the one used to construct the stencil Tg used to reconstruct a third order 



accurate divergence-free mag netic field B^?''^^) described in Sec. 13.1.21 Stencils T^^) and T^^) 
are constructed similarly. 

We next reconstruct preliminarily polynomial Pq^^ by using stencil T^^^ by solving the 
following constrained least square problem 

{ J^^Pi''\x,y)dxdy = \ICi\vi, / = 4, 5, 40, 41, 50, 51; 

subject to: (3.16) 
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Polynomials PQ^\x,y) reconstructed by using T^^^ and PQ\x,y) reconstructed by using 
T*^^) are computed similarly. 

For each P^\x^ y),m = 1, ■ ■ ■ , 7, we compute a smoothness indicator [TH] by 



1/2 



SI{Pt^)={j2[ h-\D-pt\x,y)rdxdy] . (3.17) 



This smoothness indicator is suitable for stringent shock wave interaction problems. See [18 
for discussion of other oscillation indicators. 

Weights Um from these smoothness indicators are redefined by 



/ ^ \ -4 



(e + SI{Pi 

O0m = —. ^—4, (3.18) 

Elib^m{e + SI{Pi'^)) 

where = 10 when m = 1 and 6^ = 1 otherwise, e = 10~^ is used to avoid division by 
zero. 

The final nonlinear WENO reconstruction polynomial PQ{x,y) is defined by 

7 

Po{x,y) = Y,^mPt\^^y) ■ (3.19) 

m=l 

This completes the reconstruction for approximating v{x,y) on JCq. In the present paper, 
this reconstruction algorithm is applied at the end of every Runge-Kutta stage, and used 
to reconstruct every component of U"^ with Vi replaced by the cell average values of corre- 
sponding component of computed by the base finite volume scheme. 



4 Numerical Test Problems 
4.1 Vortex evolution problem 

We consider a vortex evolution problem, which was initially suggest in [35] and was adapted 
to the MHD equations in [5], to assess the convergence order of the scheme. 

The problem is defined on a [—5, 5] x [—5, 5] domain with periodic boundary conditions on 
both sides. The unperturbed MHD flow is given by {p,Pgas,Ux,Uy, B^, By) = (1, 1, 1, 1,0,0). 
The ratio of specific heats is 7 = 5/3. The vortex is introduced through perturbed velocity 
and magnetic fields given by 

{Su.,Suy) = ^e'-'^'-^'\-y,x) , i5B,,SBy) = ^e'-^(^-^'\-y, x) , 

where r'^ = x"^ + y"^. The pressure determined by the dynamical balance is given by 

_ /.^(l-r^)-/i^ . 
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We use K = 1, /i = 1 in our computation. The exact solution is the initial configuration 
propagating with speed (1, 1), and is given by 

U{x,y,t) = Uoix -t,y- t). 

The computation domain is [—5, 5] x [—5, 5]. Periodic boundary condition is used at both side 
of the domain. The periodic boundary condition introduces an error of magnitude O{10~^), 
which does not affect the reported results. The typical triangle edge length, denoted by h, 
is listed in the first column of all the tables shown in this section. We show the Li and 
Loo errors and orders of accuracy of variables p and e at time T = 1.0. Table [T] shows the 
result for the second order accurate divergence-free WENO reconstruction-based scheme. 
The result for the third order accurate scheme is shown in Table [2j As can be seen, results 
in these two tables show clearly that we have achieved the expected accuracy property of 
the scheme. The absolute value of the undivided divergence of the magnetic field is about 
O{10~^^) in these simulations. 

Table 1: Numerical errors and convergence order for the second order accurate divergence- 
free WENO reconstruction-based method for solving the 2D vortex evolution problem. 



h 


Li 


order 


Loo 


order 


Li 


order 


Loo 


order 




p error 




p error 




6 error 




e error 




1/40 


2.93E-3 




7.65E-3 




1.07E-1 




3.69E-2 




1/80 


8.22E-3 


1.83 


4.61E-3 


0.73 


2.71E-2 


1.98 


1.42E-2 


1.38 


1/160 


1.38E-3 


2.57 


1.31E-3 


1.81 


4.95E-3 


2.46 


4.84E-3 


1.55 


1/320 


2.40E-4 


2.53 


3.92E-4 


1.74 


9.23E-4 


2.42 


1.71E-3 


1.49 



Table 2: Numerical errors and convergence order for the third order accurate divergence-free 
WENO reconstruction-based method for solving the 2D vortex evolution problem. 



h 


Li 

p error 


order 


Loo 

p error 


order 


Li 

£ error 


order 


Loo 

e error 


order 


1/20 


2.01E-2 




2.80E-3 




1.31E-1 




2.49E-2 




1/40 


2.31E-3 


3.12 


4.96E-4 


2.50 


1.03E-2 


3.67 


2.44E-3 


3.36 


1/80 


1.85E-4 


3.64 


4.44E-5 


3.48 


6.64E-4 


3.96 


2.03E-4 


3.59 


1/160 


2.58E-5 


2.84 


8.75E-6 


2.34 


9.26E-5 


2.84 


1.90E-5 


3.42 



4.2 Numerical dissipation and long-term decay of Alfven waves 

We consider a smooth solution problem proposed in [5j, which examines the numerical dis- 
sipation of torsional Alfven waves that are made to propagate at a small angle to the y-axis. 
We use the same ang le a = tan(l/6) = 9.462°; and the magnetic field is normalized by a 
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1 are initial values of density and 
: 0, and the unperturbed magnetic 



factor. The density po = 1, and pressure po - 
pressure respectively. The unperturbed velocity is uo 
field is i?o = 1- 

The computational domain is [— r/2, r/2] x [— r/2, r/2] with r = 6. The direction of wave 
propagation is along the unit vector n = (ra^,^^) = (;^/^f^, '^/^^) ■ The phase of the wave 
is taken to be = ^{n^x + Uyy — V^t) , where Va = Bq^/p^. The velocity is given by 
u = {uqiIx — eUy cos 0, UQUy + en^ cos 0, e sin 0) , where e = 0.2. The magnetic field is given 
by B = {Bqu^ + euyy/pocoscj), B^ny — en^y/po cos (j), — e^posin^) . 

The computational domain is [—3,3] x [—3,3]. The typical edge length of triangles is 
roughly equal to Solution of the problem is computed to a time T = 129. The maximum 
values of Uz and Bz should remain constant over time for the exact solution, but decay due to 
the numerical dissipation. Therefore this problem provides a good assessment of dissipation 
of the numerical scheme. Figure |4] shows the logarithm of the maximum of absolute values 
of Uz and Bz over time. We see clearly that the third order accurate scheme is substantially 
less dissipative than the second order accurate scheme. 





60 80 
time 



(a) 



(b) 



Figure 4: (a) Logarithm plot of the maximum of absolute value of the z-component of the 
velocity from the torsional Alfven waves propagation problem, (b) Logarithm plot of the 
maximum of the absolute value of the ^-component of the magnetic field from the torsional 
Alfven waves propagation problem. The maximum value should remain constant for the 
exact solution but decay due to the numerical dissipation. The solid line is for the third 
order accurate scheme; while the circled line is for the second order accurate scheme. 



In what follows, we test the problems with discontinuities to assess the non-oscillatory 
property of the proposed third order accurate scheme. 

4.3 Rotor problem 

This test problem is first proposed in ^ and is considered as the second rotor problem in 
[37] . The computational domain is [0, 1] x [0, 1]. 7 = 5/3. A dense rotating disk of fluid is 
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initially placed at the central area of the computational domain, while the ambient fluid is 
at rest. The initial condition is given by 

2.5 

{p,Pgas, Ux, Uy, M^, B^, By, B^) = (p(x, 0), 0.5, ^^^(x, 0), Uj^(x, 0), 0, — =, 0, 0, 0). 

v47r 

Here 

( 10, -(y-0.5)/ro, (x-0.5)/ro if r < tq 
(p(x,0),M,(x,0),M,(x,0))= <^ 1 + 9/, -{y - 0.5) f/r, (x - 0.5)//r if ro < r < n 

[ 1, 0, 0, if r > ri 

where ro = 0.1, ri = 0.115, f = {ri — r)/{ri — tq), and r = y^(x'^^TL5)M-~(y'--T)y5p. 

The solution at time t = 0.295 is computed. The typical edge length of triangles used 
to partition the domain is about A CFL number 0.4 is used for calculation. Figure 
[5] plots the numerical result of the density p, pressure Pgas, magnetic pressure {B^ + By)/2 
and Mach number. We see that there is virtually no diffusion of the loop's boundaries and 
no oscillations in the magnetic pressure within the loop's interior. The pressure is positive 
throughout the computational domain. The degradation in the density variable that was 
previously reported in [27] is not seen in our simulation. 



4.4 Blast wave problem 

This test problem is taken from [3]. It was about a spherical strong fast magneto-sonic 
shock propagates through a low-/3 (/3 = 0.000251) ambient plasma. We use it to show the 
advantages of the divergence-free reconstruction. The setup of the problem is as follows: on a 
computational domain [0, 1] x [0, 1], p = 1, u = 0, B^ = IOO/a/^tt, By = B^ = 0, Pgas = 1000 
within a circle centered at (0.5,0.5) of radius R = 0.1 and Pgas = 0.1 elsewhere. The final 
simulation time t = 0.01. The typical edge length of triangles used to partition the domain is 
about 0.0075. This is a stringent test problem [3]. The pressure is several orders of magnitude 
smaller than the magnetic energy. A small discretization error in the total energy can produce 
negative pressure near the shock front, as observed by others [23 121] • We used the negative 
pressure fix Strategy 1 in [2j to treat this. Briefly, in addition to evolve conservative variables 



in Eq. 2A_, we also update the entropy density on each cell in every numerical time step. 
If after reconstructing the magnetic field on a cell, the pressure computed from cell average 
values becomes negative, we derive the updated pressure from entropy and use that to form 
a new total energy density which corresponds to a positive pressure. We next use the new 
total energy density to reconstruct a polynomial approximation to the energy function; while 
density and momentum are reconstructed by using the average values computed by the base 
finite volume scheme respectively. We note that this treatment violates conservation of total 
energy locally. However, we only violate conservation in local regions by an amount that 
is smaller than the discretization accuracy. And we obtain a numerically consistent and 
positive pressure which is important for the physics of the problem. 

Figure [6] plots the numerical result of the density p, pressure Pgas, magnetic pressure 
+ By)/2 and magnitude of the velocity ^Ju1 + Uy. Owing to the large pressure placed at 
the center of the domain at the start of calculation, a strong blast wave propagates outwards, 
leaving a low density region in the center of the computational domain. We see that there 
is only minor oscillations in the density plot. Other fields are resolved nicely. 
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(c) (d) 

Figure 5: solution of the rotor problem at time t = 0.295. Thirty equally spaced contours 
are shown in each plot, (a) Density p; (b) Pressure Pgas', (c) magnetic pressure {B^ + By)/2; 
(d) Mach number. 

4.5 Orsag-Tang problem 

Here we simulate the Orszag-Tang vortex problem [29j. The initial conditions are Ux = 
-sm{y) Uy = sm{x), Ex = -sm{y), By = sin(2x), p = 7^ Pgas = J, Uz = = 0. The 
computational domain is a square [0, 2n] x [0, 2tt] with periodic boundary conditions along 
both boundaries. 7 = 5/3. The final output time t = tt. The typical edge length of triangles 
used to partition the domain is about Starting from a smooth initial condition, the fiow 
becomes very complex as expected from a transition towards turbulence gradually. Figure 
[7| shows the development of density p in the Orszag-Tang vortex problem. Also we report 
that the density and pressure have remained positive. No positivity fix was needed for this 
problem. 
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Figure 6: solution of the blast wave problem at time t = 0.01. Forty equally spaced 
contours are shown in each plot, (a) Density p; (b) Pressure Pgas', (c) magnetic pressure 
(Bl + B^)/2; (d) Magnitude of the velocity ^J^i^^^y 



5 Concluding Remarks 

In this paper we introduced a divergence-free WENO reconstruction-based finite volume 
method for solving the ideal MHD equations on two-dimensional triangular grids. The pro- 
posed method is based on the CT framework and achieves exactly divergence-free magnetic 
field. Numerical tests show that the proposed schemes have achieved the desired order of 
accuracy and the third order accurate scheme has been shown to perform very well for 
shock wave problems. While this paper only implements the second order accurate and the 
third order accurate schemes, the proposed method in principle can be generalized to three 
dimensions and to general meshes. 
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(a) (b) 




123456 123456 



(c) (d) 

Figure 7: Orsag-Tang problem. Evolution of p over time. Top left: t = 0.5; top right: 
t = 1.0; bottom left: t = 2.0; bottom right: t = 3.14. 15 equally spaced contours are used. 
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